Coagulation by Random Velocity Fields as a Kramers Problem 
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We analyse the motion of a system of particles suspended in a fluid which has a random velocity 
field. There are coagulating and non-coagulating phases. We show that the phase transition is 
related to a Kramers problem, and use this to determine the phase diagram, as a function of the 
dimensionless inertia of the particles, e, and a measure of the relative intensities of potential and 
solenoidal components of the velocity field, F. We find that the phase line is described by a function 
which is non-analytic at e = 0, and which is related to escape over a barrier in the Kramers problem. 
We discuss the physical realisations of this phase transition. 



Deutsch [lj introduced and investigated a model which 
can describe the motion of particles suspended in a ran- 
domly moving fluid. He showed that the one-dimensional 
model exhibits a phase transition between coagulating 
and non-coagulating phases as the effect of inertia of the 
particles is increased. We recently solved Deutsch's one- 
dimensional model exactly ('path-coalescence model' 0). 

This letter discusses the phase diagram for the path- 
coalescence model in higher dimensions, which is the 
most relevant case for physical applications. In the limit 
where inertial effects are negligible, the suspended parti- 
cles are advected by the flow: the theory of this limiting 
case is described in |3| ■ The case where the inertia of the 
particles play an important role is much less thoroughly 
understood. In Rcf. 4] a model is described in which 
inertial effects are treated by assuming that particles are 
advected in a perturbed velocity field. Ref. ja] gives nu- 
merical results on particle aggregation in a closely related 
model. The aggregation of buoyant particles on the sur- 
face of a turbulent liquid has also been studied Q . Here 
we give a treatment of the phase transition in higher di- 
mensions, using an exact mapping to a Kramers problem 
(the escape of particles from an attractor by diffusion). 

In the following we discuss the phase transition in two 
dimensions (Fig. The three-dimensional case is con- 
siderably harder to analyse, but the results are surpris- 
ingly similar. We derive a perturbation expansion for a 
Liapunov exponent Ai determining the phase transition, 
in powers of a dimensionless measure of the inertia, e. 
Surprisingly, the perturbative result for the phase line 
turns out to be independent of e. Our numerical simula- 
tions, by contrast, imply that the phase line does depend 
upon e. We resolve this apparent inconsistency by show- 
ing that there is a contribution 8\\ to Ai which is non- 
analytic in e, characteristic of the flux over a barrier in a 
Kramers problem: <5Ai ~ exp(— <E>/e 2 ), where <£> is the ac- 
tion of a trajectory in a Hamiltonian dynamical system. 
We conclude with a discussion of physical applications. 

We consider non-interacting spherical particles of mass 
to, radius a, and r(t) denotes the position of a typi- 
cal particle. These move through a fluid with velocity 
n(r, t) having viscosity rj. To avoid complications from 



FIG. f: Evolution of uniform initial particle distribution 
under eq. Q in the unit square with periodic boundary 
conditions, for f0 4 particles; F = f/3, e — 0.33 (coalescing 
phase), at t = (a), 10 (b), and t = 275 (c); T = 1/3, 
e = 4.30 at t = (d), 10 (e), and 275 (f). The statistics 
of the force is explained in the text, we chose C(R, At) oc 
exp[-At 2 /(2r 2 ) - R 2 /{2^ 2 )] with £ = 0.1/tt and r -> 0. 



displaced- mass effects 0, we assume that the particles 
have much higher density than the surrounding fluid. 
Our results are therefore most relevant to suspensions 
in gases, and we allow u(r,t) to be a compressible flow 
(by contrast, Refs. 0, 0| treat the case relevant to sus- 
pensions in liquids). We consider the case where the 
drag force /d r on the particle is given by Stokes' law: 
fdr = 6irr]a(u — r). The equation of motion is 



-7(r — u) 



(1) 



where 7 = Qifqa/m. Rearranging JQ) gives r — p/m, 
P = ~1P + f( r > t) where / = 777114 will be modelled by 
a random field. Linearising to obtain an equation for the 
separation (Sr, Sp) of two nearby trajectories gives 

Sr = 5p/m , Sp — —~/Sp + F(t)Sr (2) 

where F(t) is a 2 x 2 matrix with elements Fij(t) — 
dfi/drj(r(t), t). We write Sr = Xfig and Sp = YiXhg + 
Y2Xh s+7I /2, where fig is unit vector in direction 9. We 
expect that the scale variable X may increase or decrease, 
but that 9, Yi, and la approach a stationary distribution. 

The phase transition is determined by the behaviour 
of X: if the (maximal) Liapunov exponent Ai = 



2 



rn~ 1 d(log c X) /dt is negative, the particles coagulate [2|- 
Substituting the expressions for Sr and Sp into (J2J), and 
taking scalar products with fig and n e+T ^ 2 we obtain 

X = YiX/m , 9 = Y 2 /m (3) 

and 

Yi = (y 2 2 -y 1 2 )/m- 7 y 1 + J F d (t) 

y 2 = -2yiy 2 /m - 7 y 2 + F (t) (4) 

where F d (f) = h e .F(t)n e and F Q (t) = h g+7T / 2 .F(t)n e . 
From (|3J), the distribution of 9 becomes uniform on [0, 2ir] 
at large t, and the Liapunov exponent is 

Ax = (Y^/m . (5) 

The statistics of the random force is assumed to be ro- 
tationally invariant, so that the statistics of F& and F 
are independent of 9. For 9 = 0, we have = Fn 
and F Q = F21, so we obtain the statistics of Fd, F a from 
those of Fix, F21. Eqs. J1J are thus independent of 0. 
The Liapunov exponent is therefore determined by solv- 
ing a pair of coupled stochastic differential equations for 
(Yi,Y 2 ). 

To make further progress we must specify the statisti- 
cal properties of the isotropic and homogeneous random 
field /. We consider the case where (/) = 0. Without 
loss of generality we can write 

f(r,t) = V^(r,i)+VAn 3 ^,t) (6) 

where is a unit vector pointing out of the plane. 
The components of / arising from <p and ip are termed 
the potential and solenoidal components, respectively. 
In the following we assume that the correlation func- 
tion C(R, At) = (<fi(r,t)<fi(r',t')) is an even function of 
At = t — t' with support r (the correlation time), and has 
support £ (the correlation length) in R = \r — r'\. We 
also assume, for simplicity, that the fields <j> and ip are not 
correlated with each other, and have the same correla- 
tion functions, apart from a scale factor a 2 = (ip 2 ) / {(f> 2 ) . 
These assumptions are easily relaxed. 

If the correlation time r is short compared to other 
relevant time scales we can write (@J as a pair of coupled 
Langevin equations. We scale these into a dimensionlcss 
form, in which the diffusion matrix is the unit matrix 

dxi = \ — Xi + e(L2; 2 — x\)\ dt' + dwi , 

dx 2 — [ — x 2 — 2exix 2 ] dt' + dw 2 ■ (7) 

Here x % = ^JdIY, (for i = 1,2), t' = jt, (dm,) = 0, 
(dwidwj) — 2Sijdt', and 

D t = \( dt (Fa^FniO)) . (8) 

The two parameters T and e are defined as follows: 

T = D 2 /D l = (l + 3a 2 )/(3 + a 2 ) (9) 



characterises the ratio of the solenoidal and potential field 
amplitudes, | < T < 3, with T = ^ for purely potential 
fields, r = 3 for pure solenoidal force fields, and L = 1 
for equal field intensities. The second parameter 

e = D\ /2 /{ mi 3 ' 2 ) = (mC 1 ) 1 / 2 /(67rr/a) 3/2 (10) 

is a measure of the importance of inertial effects in the 
equation of motion. 

The Langevin equations Q are equivalent to a Fokker- 
Planck equation [8| 

dP 

— =TP= V.[-VP + VP] (11) 

where the advection field V has components V\ — —x\ + 
e{Tx 2 — x\) and V 2 = —x 2 — 2ex\x 2 . We write T = + 
eTi and seek a steady-state solution P satisfying TP = 0. 
The Liapunov exponent Ai is determined by computing 
(xi) with P(x), and using JSJl to obtain Ai = je(xi). 

We start by considering a perturbative approach. In 
the limit e — > the steady-state solution is 

P Q (x) = Aexp[-±(x 2 + x 2 )} = Acxp[-$ ( a; )] . (12) 

It is convenient to make a transformation of the Fokker- 
Planck operator T to a new operator TC of the form 
H = exp($ /2).Fexp(-$o/2). We write H = H + eHi 
and obtain, for 7Yo, 

Ho - (d 2 Xl - \x\) + (d 2 X2 - \x 2 ) + 1 (13) 

where (d 2 . — xf/A) are one-dimensional quantum- 
mechanical harmonic-oscillator Hamiltonians for coordi- 
nates Xi. The eigenvalues of Ho are —(k + I), where 
k and I are non-negative integers. We use a variant of 
Dirac notation denoting the eigenstates of Ho by \4>ki), 
with (x\4>oo) = exp[-$ o (a0/2], and P(x) = (x\P). In- 
troducing annihilation and creation operators, fij and dj, 
with the usual properties, such that Xi = hi + a] we find 

Hi = a\{Yxl-x\)-2a\x2Xi. (14) 

We seek a null eigenvector of H, satisfying H\Q) = 0, in 
the form of a perturbation series, \Q) — \Qo) + e\Qi) + 
e 2 IQ2 ) + -• ■ where \Qo) — \<t>oo)- The general term satisfies 
the recursion Ho\Q n +i) +Hi\Q n ) = 0. In terms of \Q), 
the expectation of Xi is given by 

(xi) = (<j>oo\ai\Q)/(<t> 00 \Q). (15) 

We find that (0oo|Q) = 1 in all orders of perturbation 
theory. Evaluating (</>oo|<2i|<3), we find that all of the 
even orders vanish, and that in the odd orders the square 
roots arising from the action of annihilation and creation 
operators cancel, so that the perturbation series for (xi) 
has integer coefficients. We obtain 

Ax=- 7 $>[(l-r)e 2 ] fc (16) 
fc>i 
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with ci = 1, C2 = 5, C3 = 60, C4 = 1105, C5 = 27120, c$ = 
828250, . . .. This predicts that the Liapunov exponent is 
positive (non-coagulating phase) for T > 1 and negative 
for r < 1, independent of e. Fig. |2t> shows that this 
surprising prediction is indeed false. 

It is possible that the Liapunov exponent may have 
a component non-analytic at e = 0, not captured by 
perturbation theory. The series IjlfiJ) is asymptotic and 
should be truncated at the smallest term, with index k*. 
The coefficients satisfy the recursion c k = (6k— 8)cfc_i + 
Yli=i c i°k-i with Co = —1/2, and we find 

c k ~ 3 fc \/2fc! (17) 

for k -> 00, so that k* ~ 1 + Int{l/[6(1 - r)e 2 ]}. We 
apply a principle expounded by Dingle 9] , which suggests 
that the error of an asymptotic series is comparable to 
its smallest term. This approach indicates that in the 
limit e — > the non-analytic term is of the form 8X\ ~ 
C exp{— 1/[6(1— r)e 2 ]} (where C might have a power-law 
dependence on e) . This still incorrectly predicts that the 
non-analytic contribution vanishes at Y = 1, so that the 
phase line is independent of e. 

We have therefore used an alternative approach: we 
write P as P = Aexp(— $) so that $ must satisfy 

V 2 $- V.V - V.V$ - (V$) 2 = . (18) 

The deterministic advection velocity field V contains 
terms which are quadratic in x. This suggests that $ 
is bounded by a cubic in s. In (JTHJ, the first two terms 
are expected to be bounded by multiples of \x\, whereas 
the remaining terms are expected to be bounded by a cu- 
bic function of \x\ . We expect that close to the origin, the 
solution is well approximated by $ ~ $ = \( x \ + #!)> 
and that far from the origin $ is well approximated by 
the solution of the equation containing only the leading- 
order terms, i.e. V.V$ + (V$) 2 = 0. This has the 
form of a Hamilton- Jacobi equation H(x, V3>) = E, 
where in our case E = and where the Hamiltonian 
is 10] H(x,p) — V(x).p + p 2 . The Hamilton- Jacobi 
equation is solved by integrating Hamilton's equations, 
then $ is obtained by integration along the trajecto- 
ries: $(x) = J dt x.p. We start classical trajectories 
infinitesimally close to the origin with initial conditions 
p = V$ = x, integrate them numerically, investigate the 
form of the trajectories, and determine the action func- 
tion &(x). We are especially interested in singularities 
of <E>, because these are expected to lead to non-analytic 
behaviour of the Liapunov exponent. 

On the ctq-axis, for P2 = 0, the vectors x, V and x 
are parallel, and p = — V. This trajectory crosses a 
singularity at the point x* = (— l/e,0), at which both 
V and p vanish, with action $(— 1/e, 0) = l/(6e 2 ). For 
< r < 2 our numerical experiments did not identify 
any other singular point with a smaller value of Con- 
sider the physical significance of this singular point in 
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FIG. 2: a A trajectory for the Langevin equations, J7J, for 
e = 0.5, r = 1, and < t' < 10 2 . The trajectory x(t) 
is most frequently located close to the origin, but occasion- 
ally it passes the unstable fixed point at ( — 1/e, 0), making 
a large excursion before returning to the origin. The non- 
analytic contribution to the Liapunov exponent is associated 
with these large excursions, b Liapunov exponent as a func- 
tion of e, for F = |(°), 1(0), 3(D): Monte-Carlo simulations of 
J2j (lines) and integration of Q (symbols), c Comparison of 
Monte-Carlo simulation of 10 (symbols) and eq. (1191 (lines), 
for T = 0.85(o), 1(D), 1.15(0) (with X = 0.34,0.28,0.215). d 
Phase diagram in two dimensions: Monte-Carlo simulations 
(o) of Q, the solid line is the function 1— V = exp[— l/(6e 2 )]. 

terms of the dynamics of a fictitious particle with po- 
sition x(t) described by the Langevin equation Q. A 
typical trajectory is plotted in Fig. The advection 
field V has an unstable fixed point at x* . To the right of 
the singularity, the particle is advected back towards the 
origin, and the probability density decreases rapidly as 
the fixed point is approached from that direction. To the 
left of the singularity, the particle is advected away, fol- 
lowing the advection field V: initially it moves to the left, 
returning to large positive x\ in a large circuit. The sin- 
gularity is therefore associated with diffusive escape from 
the attractor of V at x = 0, in which the escaping par- 
ticle is initially advected away, but returns along paths 
close to the positive x\ axis. This leads to the hypothesis 
that there should be a contribution to Ai proportional to 
exp[— $(a:*)] = exp[— l/(6e 2 )], where the prefactor may 
have an algebraic dependence on e. 

According to our numerical results, the dominant cor- 
rection to the perturbation series in the limit e — * docs 
indeed arise in this way: we find 

k* 

Ax/ 7 - X(r) exp[-l/(66 2 )] - £ c k [(l-Y)e 2 ] k . (19) 

fc=i 

We have not been able to determine the form of the pref- 
actor x analytically: it appears to be independent of e. 
In Fig. |2t we compare the Monte-Carlo simulations of 
{7J) with 1|19|) ■ Fig. |3l shows the phase diagram, as de- 



4 



termined using JSJ) and Monte-Carlo simulation of J7J. 
The computed curve is compared with an empirical fit, 
using the exponential function exp[— 1/ (6e 2 )] which char- 
acterises the escape process. 

The coefficients in also occur in the asymp- 

totic expansion of Ai'(z)/Ai(z) This suggests 

th at A x / 7 = -Re[Ai'0)/Ai(z) + -Jz\/(2^z) with z = 
— Le) _4 / 3 /4. This expression cannot be correct for 
|1 — r| < 1, because its leading order non- analytic cor- 
rection becomes smaller than the non-analytic term in 
(|19|l in the limit e — > 0. However, on setting T = 0,we 
find that it does agree with the exact solution of the one- 
dimensional problem obtained in 0], and our numerical 
results show excellent agreement with Monte-Carlo sim- 
ulations when r = 2 (where we believe that it could also 
be exact), and when L = 3. 

In summary, we have seen that the theory of coagu- 
lation by random velocity fields bears several surprises. 
We have seen that the phase transition is determined by 
the stationary state of a Langevin process. Perturbation 
theory incorrectly predicts that the phase line is inde- 
pendent of the inertia parameter e. The asymptotics of 
the high order terms of the perturbation series again do 
not predict the correct form of the non-analytic term; we 
are not aware of any other physical example where this 
procedure fails (although a mathematical example was 
suggested in [l^l)- Finally, the path-coalescence phase 
transition is driven by a non-analytic term characteris- 
ing the flux of a barrier in a Kramers problem. 

We conclude with a number of remarks and a discus- 
sion of possible applications of the effect. First, in many 
applications the suspended particles will not have equal 
masses, and it is necessary to consider how mass disper- 
sion affects the coagulation process. We have shown, in a 
perturbative framework, that two particles with masses 
differing by 8m <C m follow trajectories with separa- 
tion Ar = 5mg(t), where (\g(t)\ 2 ) remains bounded as 
t — > oo. We infer that when the masses are unequal, 
particles condense onto fragmented line segments rather 
than isolated points. The coagulation effect is therefore 
weakened, but not destroyed (Fig.^J by mass dispersion. 
Second, the structures observed in Fig. indicate that 
the area-contraction rate is much larger than the max- 
imal Liapunov exponent Ai. We have verified this by 
computing the Liapunov spectrum A4 < A3 < A2 < Ai 
of JQ) numerically (the area-contraction rate is given by 
A2 + Ai). Similar structures were observed in computer 
simulations of incrtial particles in a chaotic flow [5j. In 
two dimensions, A2 can be found from the stationary 
state of Langevin equations similar in structure to • 

Turning to specific applications, we note that the ran- 
dom velocity field u(r,t) could arise from random sound 
waves, turbulent flow, or other random disturbances. We 
do not know how to estimate L reliably for turbulent 



flow in gases. In the case of liquids where the flow is 
incompressible, our theory applies when the suspended 
particles are denser that the fluid: we predict there is no 
coagulation because L > 1. In sound waves the veloc- 
ity field is proportional to the pressure gradient, and is 
therefore pure potential: T — i, so according to Fig.|3l 
there is a coagulating phase. One possible technological 
application of the effect would be to the coagulation of 
small pollutant particles in an engine exhaust. An ultra- 
sonic noise source could increase the size of the pollutant 
particles until they are large enough to be captured effi- 
ciently by a mechanical filter. Coagulation by ultrasonic 
sources has been observed experimentally [l3|. Theo- 
retical treatments, reviewed in 0], have considered the 
case where the ultrasound has a single frequency, and the 
coagulation results from particles with differing masses 
experiencing different displacements. This letter has in- 
troduced a new mechanism for ultrasonic coagulation, 
which works even when particles have the same mass. 
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FIG. 3: Same as Fig. ^1 - c, but with masses uniformly 
distributed on an interval of half width 8m — 0.2(m). 



